function KeyRateDuration = bndkrdur(ZeroData,CouponRate,...
    Settle, Maturity,varargin)
%BNDKRDUR Key rate duration of bonds given a zero curve.
%
%   KeyRateDuration = bndkrdur(ZeroData, CouponRate,Settle, Maturity)
%
%   KeyRateDuration = bndkrdur(ZeroData, CouponRate,Settle,...
%                               Maturity,'Param1','Value1')
%
%   Optional Inputs: InterpMethod, ShiftValue, KeyRates, CurveCompounding,
%                    CurveBasis, Period, Basis, EndMonthRule, IssueDate,
%                    FirstCouponDate, LastCouponDate, StartDate, Face
%   
%   Description:
%       BNDKRDUR computes the key rate durations for 1 or more bonds given
%       a zero curve and a set of key rates (by default, the key rates are
%       each of the zero curve rates).  For each key rate, the duration is
%       computed by shifting the zero curve up and down by a specified
%       amount (ShiftValue) at that partiular key rate, computing the 
%       present value of the bond in each case with the new zero curves,
%       and then evaluating the following:
%
%       krdur_i = (PV_down - PV_up) / (PV * ShiftValue * 2)
%
%       Note that the shift to the curve is computed by shifting the
%       particular key rate by the Shift Value, and then by interpolating
%       the values of the curve in the interval between the previous and
%       next key rates.  For the first key rate, any curve values before
%       the date are equal to the Shift Value; likewise, for the last key
%       rate, any curve values after the date are equal to the Shift Value.
%
%   Inputs:
%       ZeroData - Zero Curve represented as an numRates X 2 matrix where
%                 the first column is Zero Dates and the second column is
%                 accompanying Zero Rates
%       CouponRate - numBonds X 1 vector of Coupon rates in decimal form.
%       Settle - Scalar MATLAB Date Number for the Settlement Date for all
%                the bonds and the zero data.  Note that this must be the
%                same settlement date for all the bonds and the zero curve.
%       Maturity - numBonds X 1 vector of Maturity dates.
%
%   Optional Inputs:
%       InterpMethod - Interpolation method used to obtain points from the 
%                     zero curve -- 'linear' (default), 'cubic', 'pchip'
%       ShiftValue - Value that zero curve is shifted up and down
%                    to compute duration, default is .01 (100 basis points)
%       KeyRates - Rates to perform the duration calculation for, specified
%                  as a time to maturity -- by default this is set to be
%                  each of the zero dates.
%       CurveCompounding - Compounding frequency of the curve -- default is
%                          semi-annual.
%       CurveBasis - Basis of the curve, where the choices are identical
%                    to Basis below -- default is 0 (actual/actual).
%
%   Additionally, the following bond properties can also be set via
%   parameter value pair:
%
%            Period - Number of coupons payments per year.
%                     Possible values include:
%                     0, 1, 2 (default), 3, 4, 6, 12
%
%             Basis - Day-count basis.
%                     Possible values include:
%                     0 - actual/actual (default)
%                     1 - 30/360 SIA
%                     2 - actual/360
%                     3 - actual/365
%                     4 - 30/360 PSA
%                     5 - 30/360 ISDA
%                     6 - 30/360 European
%                     7 - actual/365 Japanese
%                     8 - actual/actual ISMA
%                     9 - actual/360 ISMA
%                    10 - actual/365 ISMA
%                    11 - 30/360 ISMA
%                    12 - actual/365 ISDA
%                    13 - bus/252
%
%      EndMonthRule - End-of-month rule; default is 1 (in effect)
%                     0 - Rule is NOT in effect for the bond(s)
%                     1 - (default) Rule is in effect for the bond(s) (meaning
%                         that a security that pays coupon interest on the last
%                         day of the month will always make payment on the last
%                         day of the month)
%
%         IssueDate - Bond issue date.
%
%   FirstCouponDate - Irregular or normal first coupon date.
%
%    LastCouponDate - Irregular or normal last coupon date.
%
%         StartDate - Forward starting date of payments.
%
%              Face - Face value of the bond; default is 100.  Note that
%                     Face has no impact on key rate duration -- this
%                     calling sequence is preserved for consistency.
%
%   Outputs:
%       KeyRateDuration - [numBonds x numRates matrix] of the 
%                         key rate durations.
%
%   Example:
%
%   This example is based on example from Golub and Tilman:
%
%       ZeroRates = [0.0476 .0466 .0465 .0468 .0473 .0478 ...
%                   .0493 .0539 .0572 .0553 .0530]';
%
%       ZeroDates = daysadd('31-Dec-1998',[90 360 360*2 360*3 360*5 ...
%                   360*7 360*10 360*15 360*20 360*25 360*30],1);
%       ZeroData = [ZeroDates ZeroRates];
%       krdur = bndkrdur(ZeroData,.0525,'12/31/1998',...
%                   '11/15/2028','KeyRates',[2 5 10 30])
%
%   See also BNDDURY, BNDCONVY, BNDCONVP.

%   References:
%   [1] Tuckman, B. Fixed Income Securities. Hoboken, NJ: John Wiley &
%   Sons, Inc., 2002.
%   [2] Golub, B.W. and Tilman, L.M. Risk Management: Approaches for Fixed
%   Income Markets. Hoboken, NJ: John Wiley & Sons, Inc., 2000.

%   Copyright 1995-2008 The MathWorks, Inc.
%   $Revision: 1.1.6.4 $   $Date: 2009/05/07 18:23:20 $

% Parse input arguments
p = inputParser;

p.addParamValue('interpmethod','linear',@(x) any(strcmpi(x,{'linear','cubic','pchip'})));
p.addParamValue('shiftvalue',1e-2);
p.addParamValue('keyrates',[]);

% Parse the bond inputs with no checking, do that in INSTARGBOND
p.addParamValue('period',[]);
p.addParamValue('basis',[]);
p.addParamValue('endmonthrule',[]);
p.addParamValue('issuedate',[]);
p.addParamValue('firstcoupondate',[]);
p.addParamValue('lastcoupondate',[]);
p.addParamValue('startdate',[]);
p.addParamValue('face',[]);

p.addParamValue('curvecompounding',2,@(x) ismember(x,[-1 1 2 3 4 6 12]));
p.addParamValue('curvebasis',0,@(x) isvalidbasis(x));

if ~isnumeric(ZeroData) || (size(ZeroData,2) ~= 2)
    error('finance:bndkrdur:zeroDataNumeric',...
        'Zero Data must be a numeric matrix with 2 columns.');
end

ZeroDates = ZeroData(:,1);
ZeroRates = ZeroData(:,2);

% Check to make sure that Settle is a scalar
Settle = datenum(Settle);
if length(Settle) > 1 || ~isnumeric(Settle)
    error('finance:bndkrdur:settleError',...
        'Settle must be a scalar MATLAB date number')
end
 
CurveSettle = Settle;

try
    p.parse(varargin{:});
catch ME
    newME = MException('finance:bndkrdur:optionalInputError',...
        'Error in optional parameter value inputs');
    newME = addCause(newME,ME);
    throw(newME)
end

[CouponRate, Settle, Maturity, Period, Basis, EndMonthRule, IssueDate, ...
    FirstCouponDate, LastCouponDate, StartDate, Face] = ...
    instargbond(CouponRate,Settle,Maturity,...
    p.Results.period, p.Results.basis, p.Results.endmonthrule, ...
    p.Results.issuedate,p.Results.firstcoupondate,  ...
    p.Results.lastcoupondate, p.Results.startdate,  p.Results.face);

InterpMethod = lower(p.Results.interpmethod);
ShiftValue = p.Results.shiftvalue;
CurveCompounding = p.Results.curvecompounding;
CurveBasis = p.Results.curvebasis;
KeyRates = p.Results.keyrates;

% Get the cash flows
[CFlowAmounts, CFlowDates] = cfamounts(CouponRate,...
    Settle, Maturity, Period, Basis, EndMonthRule, IssueDate,...
    FirstCouponDate, LastCouponDate,StartDate, Face);

nBonds = length(CouponRate);

CFTimes = zeros(size(CFlowDates));
CFRates = zeros(size(CFlowDates));
DF = zeros(size(CFlowDates));

CurveDF = zero2disc(ZeroRates, ZeroDates, CurveSettle, CurveCompounding, CurveBasis);

% This loop is needed, because each bond's compounding and basis can be
% different
for bondidx=1:nBonds
    if isisma(Basis(bondidx))
        BondCompounding = 1;
    else
        BondCompounding = 2;
    end
    ZeroRatesNew = disc2zero(CurveDF,ZeroDates,CurveSettle,BondCompounding,...
                                                        Basis(bondidx));
    
    ZeroTimes = yearfrac(CurveSettle,ZeroDates,Basis(bondidx));
    CFTimes(bondidx,:) = yearfrac(CurveSettle,CFlowDates(bondidx,:),...
                                            Basis(bondidx));
    
    CFRates(bondidx,:) = interp1(ZeroTimes,ZeroRatesNew,...
                        CFTimes(bondidx,:),InterpMethod,'extrap');
    
    notnanidx = ~isnan(CFlowDates(bondidx,:));
    DF(bondidx,notnanidx) = zero2disc(CFRates(bondidx,notnanidx),...
        CFlowDates(bondidx,notnanidx),CurveSettle,BondCompounding,Basis(bondidx));
end

% Compute current prices
PV = nansum(CFlowAmounts.*DF,2);

if isempty(KeyRates)
    KeyRates = yearfrac(CurveSettle,ZeroDates,CurveBasis)';
end

nRates = length(KeyRates);

KeyRateDuration = zeros(nBonds,nRates);

% Compute key rate durations
for rateidx=1:nRates
    
    KeyRateShift = zeros(1,nRates);
    KeyRateShift(rateidx) = ShiftValue;
    ShiftRates = interp1(KeyRates,KeyRateShift,CFTimes,InterpMethod,0);
    
    if rateidx==1
        earlyidx = CFTimes < KeyRates(rateidx);
        ShiftRates(earlyidx) = ShiftValue;
    elseif rateidx==nRates
        lateidx = CFTimes > KeyRates(rateidx);
        ShiftRates(lateidx) = ShiftValue;
    end
    
    UpRates = CFRates + ShiftRates;
    DownRates = CFRates - ShiftRates;
    
    % Loop through each bond to handle different basis and compounding
    DF_up = zeros(size(DF));
    DF_down = zeros(size(DF));
    for bondidx=1:nBonds
        
        if isisma(Basis(bondidx))
            BondCompounding = 1;
        else
            BondCompounding = 2;
        end
        
        notnanidx = ~isnan(CFlowDates(bondidx,:));
        
        DF_up(bondidx,notnanidx) = zero2disc(UpRates(bondidx,notnanidx),...
            CFlowDates(bondidx,notnanidx),CurveSettle,BondCompounding,Basis(bondidx));
        DF_down(bondidx,notnanidx) = zero2disc(DownRates(bondidx,notnanidx),...
            CFlowDates(bondidx,notnanidx),CurveSettle,BondCompounding,Basis(bondidx));
    end
    
    % Compute PV_up and PV_down from the cash flows and discount factors.
    PV_up = nansum(CFlowAmounts.*DF_up,2);
    PV_down = nansum(CFlowAmounts.*DF_down,2);
    
    % Compute duration
    KeyRateDuration(:,rateidx) = (PV_down - PV_up)./(2*ShiftValue*PV);
end
